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ABSTRACT 

We discuss the implementation of a new regular algorithm for simulation of the gravitational 
few-body problem. The algorithm uses components from earlier methods, including the chain 
structure, the logarithmic Hamiltonian, and the time-transformed leapfrog. This algorithmic 
regularization code, AR-CHAIN, can be used for the normal iV-body problem, as well as for 
problems with softened potentials and/or with velocity-dependent external perturbations, in- 
cluding post-Newtonian terms, which we include up to order PN2.5. Arbitrarily extreme mass 
ratios are allowed. Only linear coordinate transformations are used and thus the algorithm is 
somewhat simpler than many earlier regularized schemes. We present the results of performance 
tests which suggest that the new code is either comparable in performance or superior to the 
existing regularization schemes based on the Kustaanheimo-Stiefel (KS) transformation. This is 
true even for the two-body problem, independent of eccentricity. An important advantage of the 
new method is that, contrary to the older KS-CHAIN code, zero masses are allowed. We use our 
algorithm to integrate the orbits of the S stars around the Milky Way supermassive black hole 
for one million years, including PN2.5 terms and an intermediate-mass black hole. The three S 
stars with shortest periods are observed to escape from the system after a few hundred thousand 
years. 



Subject headings: black hole physics 
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1. Introduction 

After the introduction of electronic computers, 
those who carried out simulations of the gravita- 
tional A^-body problem soon realized that the clas- 
sical methods of numerical integration were often 
not satisfactorily accurate due simply to the strong 
l/r^ character of the gravitational force. 

Th e situation changed when lKustaanheimo fc Stiefell 
(jl965l ) published their (KS-) transformation from 
four-dimensional to three-dimensional space. The 



special case of planar system h ad been known for 
a long time ( Levi-Civita 1923) but it had turned 
out that a similar transformation in the three- 
dimensional space was not possible. The situation 
for the general A^-body problem improved only af- 
ter the KS transformation beca me well known, due 
i n larg e part to publication bv lStiefel fc Scheifele 
( 19711 ) of their text, which comprehensively dis- 
cussed the application of the KS-transformation 
to the perturbed two-bod y proble r n. F u rther 
on, lAarseth fc Za^ (|l974l) . iHeggid (|l974l) . IZard 



(|l974l ) and lMikkola &: AarsethI (|l993l ) applied the 
KS transformation to the general three-body and 
A^-body problems. 

An entirely new way of regularizin g close en- 
counters was inve nted simultaneously bylMikkola & 
|l999a,b) and bv lPreto fc Tremaind (fl999 1. This 
new method introduced the so-called logarithmic 
Hamiltonian (LogH). Together with the simple 
leapfrog algorithm, this new method gives reg- 
ular results for close encounters; in fact a cor- 
rect trajectory is obtained for the two-body prob- 
lem. Other, similar meth ods are described by 
Huang fc Leimkuhleil (|l997t ). 

The remaining problem was that none of these 
regularization methods could be easily applied to 
systems with extremely large mass ratios. In an 
attem pt to solve this problem, Mikkola fc Aarse"tlil 
(|2002[ ) introduced the time-transformed leapfrog 
(TTL). This method is in some cases mathemat- 
ically equivalent to the LogH method, but it is 
more general, and arbitrary mass ratios are al- 
lowed. The drawback of this method, however, 
is that in some cases the roundoff error affects the 
results considerably. In these new methods (LogH, 
TTL) the regularization is achieved by using the 
leapfrog, hence the name "algorithmic regulariza- 
tion." More details about the thus-far mentioned 
metho ds can be found in the book by lAarsethl 
(|2n03h . 

Since the leapfrog alone is rarely accurate 
enough, one must supplemen t the metho d with 
the extrapolation algorithm (Grag^ jl96i 119651: 
Bulirsch fc Stoeilll966t IPress et al.lll 986'). or with 
higher-order leapfrogs (jYoshid J|l990), in order 
to get highly accurate results. Efhciency of the 
extrapolation procedure requires that the basic 
algorithm (leapfrog in algorithmic regularization, 
modified midpoint method in the KS-regularized 
codes) have a certain symmetry. In the case of 
the leapfrog this means time reversibility. If the 
system has velocity-dependent forces, such as rela - 
tivistic post-Newtonian (PN) terms (|Soffellll989l ). 
then the required symmetry is more difficult to ob- 
tain. One w ay to cure this prob l em w as recently 
obtained by iMikkola fc MerrittI (|2006h . who for- 
mally doubled the dimensionality of the parame- 
ter space and constructed a generalized midpoint 
method (GAR). This algorithm allows the use of 
algorithmic regularization and velocity-dependent 
perturbations. However, these authors discussed 



the new algorithm essentially only for the case of 
the perturbed two-body problem. 

In this paper we discuss the details of our most 

rece nt implementation of the algorithmic regular- 

Tanikaiwalion method. This uses the chain structure, the 
same structure as in the KS-CHAIN algorithm of 



Mikkola fc AarsethI (|l993[ ). That device, together 



with a new time-transformation function, signifi- 
cantly reduces the roundoff problems and makes 
the new code a good alternative for simulations 
of strongly-interacting few-body systems. In the 
final section, we present some applications of the 
chain routine to the problem of relativistic orbits 
around the supermassive black hole at the Galac- 
tic center. 

2. Notation 

In this paper, we use the following basic sym- 
bols: 

G — 1 gravitational constant 

t time 

nik mass of body k 

Pfc momentum 

Tfc position 

Vfc = r'fc velocity 

T = Y,k=i Pfe/2wfe kinetic energy 

U = J2o<i<j<N rriimj/lri - rj\ potential energy 

E = T — U total energy 

B = —E = U — T binding energy 

3. Regular Algorithms 

In the few-body problem, close approaches of 
two bodies are common, and these events re- 
quire highly accurate computations since the en- 
ergies involved are large. Thus a method that 
can accurately advance the motions of few-body 
systems must be accurate and efficient for the 
two-body problem. In addition to the classical 
KS-transformation there are some new algorithms 
that satisfy this requirement. 

An additional, and most important, problem 
is the roundoff error. This becomes a major prob- 
lem if e.g. center-of-mass coordinates are used and 
there are close binaries and/or close encounters in 
the simulated system. The cure for this problem 
is the use of the chain structur e, originally applied 
with the KS-transformation ( Mikkola fc AarsethI 
19931 ) in order to KS-regularize all the short(est) 
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distances. Later it became clear that the chain 
structure was also beneficial in reducing signifi- 
cantly the roundoff error. 

In this section, we first review the ingredients 
we need in the construction of the final A'^-body 
algorithm. These basic algorithms produce reg- 
ular results in close approaches and their results 
can be improved to hig h precision using e.g. the 
Bulirsch fc Stoei (Il966h extrapolation method. 



3.1. LogH, TTL and GAR 

3.1.1. LogH 

Recently, iMikkola fc Tanikawal (|l999a[ l and 
IPreto fc Trem ainc (1999) pointed out that the 
logarithmic Hamiltonian in extended phase space: 



A = ln(T + S) -ln([/). 



(1) 



where B (the binding energy) is the momentum of 
time, gives the equations of motion in the form 



t' = 




(2) 


r', = 


dT 


(3) 


B' = 




(4) 


Pfe = 




(5) 



Here we include equation @ although this is 
needed only if there is a time-dependent poten- 
tial to be added to the A^-body potential U . Since 
the right hand sides of these equations do not de- 
pend on the left hand side variables, a leapfrog 
algorithm is possible. That may be symbolized as 



X(/i/2)V(/i)X(/i)..V(/i)X(V2), 



(6) 



where X(s) means solution for the coordinate 
equations ([2|), ^ with constant B and pfc over 
an integration step of length^ s: 



St::^ s/{T + B); t^t + 6t: rk ^ Vk + St 



dT 
dpk' 



(7) 



Correspondingly V(s) signifies the operation 



5t = s/U; B ^ B + 6t 



dU 



Pk Pfc + St 



dU 
dvi, ' 



(8) 



which solves equations (g]) and (O for constant t 
and Tk. 



For the case of two bodies only, this algorithm 
produces correct trajectories, with only an OQi^) 
phase error. This is true even for the collision 
orbit and the energy conservation is, in this case, 
typically in the machine precision level. 

3.1.2. TTL 

The second important ingredient i n our algo- 

rithm is the time-transformed leapfrog ( Mikkola fc Aarse"thl 
20021 ). The basic idea is to introduce a coordinate- 
dependent time transformation function $!(.., Yk, ..) 
and a new variable oj which is supposed to have 
the same numerical value as $7, but that value is 
obtained from the differential equation 



^ dvk 

k 



(9) 



The equations of motion can be written as 

t' = 1/lo (10) 

^'k = Vfc/^ (11) 



/o 

^ OVk 



(12) 
(13) 



where is the acceleration of the /c'th particle. 
The structure of these equations allows the con- 
struction of a leapfrog algorithm: 

X(s) 

5t = s/lo\ t^t + St; Tk^Tk + St iU) 



V(s) 

6t = s/ri; Vfc 



Vfc 



StAk 



(15) 

10 + St <iM] 



where < 51 > is the average of this quantity over 
the step, i.e. 



u , v^;""')/2. (17) 



Here the superscripts "old" and "new" refer to Vfc 
values before and after the velocity advancement 
in the operation (|16p. 

3.1.3. GAR 

Here we concisely review the G AR method. 
Following iMikkola fc MerrittI (|2006l ). we consider 
a differential equation 



z = f (z), 



(18) 
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and an approximation for its solution (over a short 
step = h) written as 



z{h) w z(0) +d(z(0),/i). 



(19) 



Here the increment d{zQ,h) can be any approx- 
imation suitable for the particular equation in 
question. If one writes, instead of p8)) . the two 
equations 

x-f(y); y-f(x), (20) 

and solves this pair with the initial values x(0) = 
y(0) = 2(0), the solution obviously is x(<) = 
y{t) = z(t). Using the pair of equations one can 
write a leapfrog as 

XI = xo + ■^f(yo); yi = yi + hi{xi); xi = xi + 

which is actually nothing but the well known mod- 
ified midpoint method. In the above one can split 
the advancement of y in two operations to get 

XI = xo + ^f(yo); yi =yo + ^f(xi)(21) 
yi = yi + Sf(xi); xi=xi + J/(yi)(22) 
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2 2' 



yi = yi 



d(xi, -); xi = XI - d(yi,-i 



encounters if the coordinates of the approaching 
bodies are measured from a distant origin. This 
problem is significantly reduced by utilizing the 
chain structure. This was originally used in a 
KS-r egularization algorithm (jMikkola &: Aarseth 



19931 ) in order to get all the short distances reg- 



ularized by the KS-transformation, but here the 
sole purpose is the roundoff reduction, for which 
the device has proved itself. 

4.1. Basic Formulation 

Let r = -i J2k=i "^fc'^fc be the kinetic energy, 
and U — ruimj /rijt\ie potential such 

that the total energy is E = T — U; the binding 
^ energy is B = U — T. 

One forms a chain of particles such that 
the shortest relative vec tors are in the chain 
( Mikkola fc Aarseth 1993f) . We stress again that 
the main purpose of using the chain structure in 
this method is to reduce the (often significant) 
effects of roundoff error. 

Let us collect the chain coordinates Xfe — r,;^ — 
in the vector X = (Xi, X2, .., Xjv-i) and let 



and this can be readily generalized using the more 
general increment d(z, h). This results in the gen- 
eralized midpoint method 



XI = xn-Fd(yo,-); yi =yo-d(xi,--|^3) 



This method has the great advantage that one 
can use any special approximation d and the al- 
gorithm is time reversible (albeit only in the ex- 
tended xy-space) and thus suitable for being used 
as the basic integrator in an extrapolation method. 
Specifically one may stress that the increment d 
may be computed using LogH or TTL which pro- 
duce good approximations in case of close en- 
counters. In addition, appearance of velocity- 
de pendent forces is here not problematic, as shown 
bv lMikkola fc MerrittI (|200e 



4. The AR-CHAIN Algorithm 

While in the algorithms discussed above one 
can use any coordinate system, the roundoff er- 
ror may be a major problem in the case of close 



•■3k 



Vjfc be 



the corresponding velocities = Vj 
in the vector V = (Vi, V2, .., Vat-i). Then the 
Newtonian equations of motion may be formally 
written 



X = V 
V = A(X) 



(25) 
(26) 



where A is the N-body acceleration and f is some 
external acceleration (e.g. due to other bodies). 

One may use the two equivalen t time transfor- 
mations (|Mikkola fc MerrittI [2006h 



ds = [a{T + B)+ l3uj + -f]dt = [aU + fdfl + j]dt, 

(27) 

where s is the new independent variable, a, (3 and 
7 are adjustable constants, is an optional func- 
tion of the coordinates = ri(X), while the initial 
value u!{0) = il(0) and the differential equation 

determine the values of uj (in fact, uj{t) — fl(t) 
along the exact solution) . 

The time transformation thus introduced regu- 
larizes the two-body collisions if one uses the sim- 
ple leapfrog algorithm as a basic integrator; re- 
sults from which can, and must, be improved using 
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an extrapolation meth od (e.g. iBulirsch fc Stoei 
19661 : IPress et al.lll986l ). 



It is possible to divide the equations of motion 
into two categories (where derivatives with respect 
to the new independent variable s are denoted by 
a prime). 

Coordinate equations: 



t' = l/{a{T + B)+puj + j) 
X' = t' V 

Velocity equations: 



(29) 
(30) 





= l/{aU + f3n + j) 


(31) 


v 


= ?(A + f) 


(32) 


/ 

UJ 


- i'^-V 
dX 


(33) 


B' 


^ dV 


(34) 



In these equations the right hand sides do not de- 
pend on the variables on the left. Consequently 
it is possible to construct a regular leapfrog algo- 
rithm for obtaining the solutions. The leapfrog 
results then can easily be improved with the ex- 
trapolation method. 

4.2. Details 

4-2.1. Finding and updating the chain 

First we find the shortest interparticle vector 
which is adopted as the first part of the chain. The 
chain is then augmented by adding the relative 
vector to the particle nearest to one or the other 
end of the existing chain. When all particles are 
included, they are re-numbered along the chain as 
1,2, ..N for ease of programming. 

To reduce roundoff problems, the transfor- 
mation from the old chain vectors Xfc to the 
new ones is done directly by expressing the new 
chain vectors as sums of the old ones as in 
Mikkola fc AarsethI (|l993h . 



4.2.2. Transformations 

When the particles are renamed along the chain 
as 1, 2, . . . , A'' one can evaluate 



Xfc = 



- Vfc. 



(35) 
(36) 



The center-of-mass quantities are 

M = ^mfe (37) 

fc 

Tern = ^mkTk/M (38) 

fc 

Vcm = ^mfcVfc/M. (39) 



The inverse transformation is done by simple sum- 
mation 



?i = 

vi = 

Tfc+i = Tk + Xfc 

Vfc+i = Vfc + Vfc, 



followed by reduction to the center of mass 

fern = ^ mkYk/M 
k 

Vcm = y^TOfcVfc/M 



rfe 

Vfc 



k 

r fc ■ 

Vfc 



- Vr 



(40) 
(41) 
(42) 
(43) 



(44) 

(45) 

(46) 
(47) 



Note that it is not always necessary to reduce the 
coordinates to the center-of-mass system since ac- 
celerations only depend on the differences. 

4.2.3. Equations of motion and the leapfrog 
One writes the equations of motion as 

Xfc = Vfc (48) 
Vfc = Ffc.+1-Ffc+ffc.+i-ffc, (49) 

where ffc are the individual external accelerations 
and the iV-body accelerations Ffc are 



(50) 



where, for j < k 

Tfc 



if fc > j -I- 2 
r,fc = { Aj-; if fc = j + 1 
Xj + Xj+i; if fc = j + 2 



(51) 



For k > j one uses Vjk — —J^kj- In the acceleration 
computation the use of Xj and Xj +Xj+i reduces 
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the roundoff effect significantly. This is one of the 
most important features of the algorithm. 

The kinetic energy is evaluated as usual 



while the potential 



(52) 



(53) 



is obtained along with the accelerations according 
to ([STjl . For the time transformation function Q. it 
seems advantageous to use 



(54) 



Here fiij are adjustable constants (see below). 

Now one is able to evaluate the two time trans- 
formation functions 

t' = l/(a(T + B) +/3t^ + 7) (55) 
t' = l/(aC/ + /31} + 7), (56) 

which are equivalent along the correct solution i.e. 
t' = i'. The evolution of with w(0) = 1^(0), is 
obtained by 

(57) 



k 



In the presence of external perturbations the bind- 
ing energy evolves according to 



B 



E 

k 



(58) 



The leapfrog for the chain vectors and Vfe can 
be written as the two mappings 

Ms): 

St = s/{a{T + B)+PLj + -f) (59) 
t = t + St (60) 
Xfc ^ Xk+dtVk (61) 

(62) 



Vis) 



St = s/{aU + pn + -f) (63) 
Vfc ^ Vk + St{Fk+i-Fk + {k+i-m) 
B ^ B + StY, i-rrik < Vfe > -ffe) (65) 



k 



^ + ^* } .T, <^k >, 

art 

k " 



(66) 



where < v/j > is the average of the initial and 
final v's (obtained from the V's according to the 
equations in section (|4.2.2p . It is necessary to eval- 
uate the individual velocities because the ex- 
pressions for B' and oj' in terms of the chain vector 
velocities are rather cumbersome. 
A leapfrog step can be written as 

X(V2)V(/i)X(V2) 

and a sequence of n steps as 

X(V2) [KZ\{^{h)l^{h))] V(/.)X(V2). 

This is useful with the extrapolation method when 
advancing the system over a total time interval of 
length = nh. 

4.3. Time Transformation Alternatives 

If one takes 



Vtj = mirrij, 



(67) 



then a — 0, (3 — 1, J — is mathematically 
equivale nt to a = 1 , /? = 7 = as was shown in 
.Mikkola fc AarsethI ([2002). However, numerically 
these are not equivalent, mainly due to roundoff 
errors in updating the value of fi, and the LogH 
alternative is numerically more stable. However, 
for proper treatment of small bodies some function 
VI is to be used. 

For increased numerical stability in the motions 
of the large bodies, and smoothing of the encoun- 
ters of small bodies, the recipe is a = 1, /3 7^ 
(but small) and 



if rrii * irij < e vn? 



= m 
= 0; otherwise 



(68) 



where = Ej<j ~ l)/2) is the 

mean mass product and e ~ 10~'^. It may be 
advisable to integrate ([57)1 for w even if /3 = 0, 
in order to force the the integrator (extrapolation 
method!) to use short steps if Co is large, thus 
giving higher precision when required. In fact, nu- 
merical experiments suggest that in most practical 
cases the parameters (a,/3, 7) — (1,0,0) give the 
best results. Exceptions are cases with extremely 
large mass ratios such that the contribution of the 
small masses in the potential are negligible (e.g. 
zero masses included). This point, however, needs 
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further investigation. One potential problem in 
the integration of the quantity w is that the incre- 
ments of it can be arbitrarily large if a collision of 
point masses occur. In this case the roundoff er- 
rors in the value of uj are significant, but do not af- 
fect much the results if the value of /3 is small. One 
should also realize that the parameters (a,/?, 7) 
can be changed during the integration (but not 
during an Bulirsch-Stoer integration step). 

4.4. Velocity-Dependent Perturbations 

For the case of velocity-dependent perturba- 
tions f = f(X, V), which occur e.g if one in- 
troduces relativistic post-Newtonian terms, re- 
lated algorithms were d iscussed in detail by 
Mikkola fc MerrittI (|2006l ) (although mostly for 
the perturbed two-body problem). Here we 
present those ideas in the short notation used 
in equation ([26ll . 



Implicit midpoint method: We have 

V'=f(X)(A(X)+f(X,V)), (69) 

which should be solved for constant t, X. 
This is often not easy, but it is possible to re- 
place the exact solution by the implicit mid- 
point method, i.e. one solves (often itera- 
tively) the increment of V from 

AV = 5t(X)(A(X) +f(X,Vo + iAV), 

(70) 

where Vq is the value of V before the up- 
date V — > V -f AV. This operation thus 
replaces the one in (l64|) . and the mean ve- 
locities (needed in updating B and uj) are 
obtained from Vq + ^AV. 

Generalized midpoint method: When the gen- 
eralized midpoint method is used, instead of 
the leapfrog, one can use a leapfrog step to 
obtain the increments d(X,V, s) and here 
one can use for the velocity the most re- 
cent available value of V. One thus simply 
evaluates 

AV-^t(X)(A(X)+f(X,Vo)). (71) 

The formulation of the method then guar- 
antees that the final approximation has the 
correct symmetry (and the correct form of 
the error expansion) so that the efficient 



Bulirsch-Stoer extrapolation method can be 
used. 

However, one question to be addressed here 
is: which of the two above methods is most 
efficient? This is problem-dependent and 
numerical experiments may be necessary to 
answer the question. In the implementa- 
tion of our code we start with the implicit 
midpoint method (which is most efficient if 
the number of iterations remains small) , and 
compute a long time average of the required 
number of iterations Q recursively as 



0.99Qo;d + -OlQn 



(72) 



This slowly "forgets" the old values and does 
not grow too quickly if there are occasional 
cases that require many iterations. However, 
when the average number of iterations ex- 
ceeds some limit, then there is time to switch 
to the GAR that does not require iterations 
(but is otherwise more expensive). 

5. Numerical Demonstrations 

In this section, we demonstrate the perfor- 
mance of the new algorithm and compare it with 
the celebrated KS- t ransfo rmed CHAIN code of 
Mikkola fc AarsethI (|l993h . We ffist considered 



a test problem consisting of one massive particle 
(the "black hole") with m — 1, and seven addi- 
tional particles ("stars") with masses 10^^ < m < 
10~^. Relativistic terms were not included. Initial 
velocities were all zero so the stars moved initially 
on nearly rectilinear orbits toward the black hole, 
but their orbits become eccentric ellipses as they 
experience perturbations from the other stars. We 
set (a,/3,7) = (1,0,0). 

We integrated the initial values given in Table 1 
including 2, 3, ..,8 bodies and carrying out integra- 
tions up to time 1000 (G — 1) using both the old 
KS-CHAIN code and the new AR-CHAIN. Two 
sets of experiments were conducted, one in which 
the Bulirsch-Stoer integrator automatically chose 
the stepsizes (in a supposedly optimal way) and 
the second with (iteration to) exact output times 
of interval At ~ 0.1. 

The results are concisely summarized in Tables 
2 and 3. There we give for both methods the exe- 
cution times and average values of the errors. The 
errors are defined in terms of the Hamiltonians 
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100 200 300 400 500 600 700 800 900 1000 1100 

time 



Fig. 1. — Evolution oilog[{T+B)/U] and 1.5*dE/L in the integration of the initial conditions in Table 1. The 
two lines are for the two different methods: the present method (AR-CHAIN) and the KS-regularized chain 
method (KS-CHAIN). In both computations the codes were allowed to choose automatically the stepsize. 



Table 1: Initial conditions for the integration of Figure 1. 



m 




X 


Y 


Z 


Vx 


Vy 


Vz 


1 






















1 X 10" 


3 


-.432498862 


-.765730892 


-.432498862 











1 X 10- 


4 


.534279612 


.288862435 


.534279612 











1 X IQ- 


5 


-.383536991 


.601722629 


-.383536991 











1 X 10" 


6 


.233942789 


-.166401737 


.233942789 











1 X 10- 


7 


.703086026 


.748854732 


.703086026 











1 X 10" 


8 


.449061307 


.186538286 


.449061307 











1 X 10" 


9 


.320791289 


-.848655159 


.320791289 












Table 2: Automatic (free) stepsize 



AR-CHAIN 


KS-CHAIN 


NB 


sec 


{\dE/U\) 


sec 


{\1.5*dE/L\) 


2 


0.28 


9.0e-13 


0.6 


4.9E-13 


3 


3.9 


1.8e-13 


2.8 


4.0E-13 


4 


8.1 


2.6e-13 


8.9 


4.1E-13 


5 


27.2 


5.0e-13 


67.9 


5.8E-13 


6 


43.6 


5.5e-13 


61.8 


1.2E-12 


7 


59.5 


l.Oe-12 


638.7 


2.9E-12 


8 


60.5 


3.0e-12 


2730.0 


7.8E-11 
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Table 3: dt=0.1 (iteration to exact output-time) 



AR-CHAIN 


KS-CHAIN 


NB 


sec 


{\dE/U\) 


sec 


{\1.5*dE/L\) 


2 


1.3 


7.5E-14 


2.2 


3.0E-12 


3 


5.9 


6.3E-13 


6.1 


1.9E-12 


4 


12.4 


3.8E-13 


14.0 


2.4E-12 


5 


32.9 


7.6E-13 


49.6 


2.2E-12 


6 


38.8 


7.3E-13 


80.5 


3.1E-12 


7 


56.6 


7.4E-13 


86.8 


2.9E-12 


8 


76.5 


9.4E-13 


12230.4 


1.3E-10 



associated with the methods: for AR-CHAIN the 
average was taken over | log((r — E)/U\ w 5E/U 
and for KS-CHAIN 1.5|(r-J7-i;)/L| « 1.56E/L. 
Here L = T+U=the Lagrangian and the factor 1.5 
is included because in virial equilibrium L — 1.5U. 
On the other hand the system is chaotic and the 
solutions are, after some time, often different (as 
with any method) and so the results must be con- 
sidered as giving just a general view of the perfor- 
mance. 

One sees that the accuracies are comparable, 
with some advantage in favour of the new method. 
Especially we note that the new method is faster 
for the case of iV = 2, i.e. a binary with eccen- 
tricity e = 1 while the precision is essentially the 
same. Thus one can conclude that this method 
works equally well, and perhaps better, than the 
celebrated KS-transformation. This is a conse- 
quence of the fact that for the two-body system 
the logarithmic Hamiltonian -I- leapfrog produces 
an exact trajectory with only a small error in the 
time, even for the collision orbit with eccentricity 
e = 1. 

The accuracy of the angular momentum is typ- 
ically similar to that of the energy. This is largely 
due to the fact that the basic algorithm (leapfrog) 
used in the method conserves angular momentum 
exactly. 

However, the execution times for cases in which 
smaller and smaller bodies are included differ con- 
siderably in favour of the new AR-CHAIN code. 
This is not very surprising, since zero masses are 
a singularity for the KS-CHAIN, but AR-CHAIN 
can handle zero masses too. 

Figure [1] compares the energy conservation for 
the two different methods in the 8-particle inte- 
gration from the initial conditions in Table 1. The 




t (yr) 

Fig. 2. — Periastron advancement of a star around 
the Milky Way supermassive black hole. The semi- 
major axis is 0.01 pc and three different values of 
the eccentricity were tried, e = (0.9,0.98,0.99). 
Dots show the major axis orientation each 1000 
orbits while solid (red) lines show the PN2.5 pre- 
diction, equation ([73|) . 



9 



I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I I 




Fig. 3. — Integration of a star on an eccentric orbit 
around an IMBH/SMBH binary at the Galactic 
center. The initial orbit elements a re sim ilar to 

2005h . 



those of the star SO-16 (Ghez et al 



Top 

panel shows the semi-major axis of the star's orbit 
with respect to the massive binary, lower panel 
shows the advancement of the periastron. The red 
(large) filled dots in this panel are the prediction of 
equation ([73]) assuming that a and e remain fixed 
at their initial values. The blue (small) dots show 
the predicted advancement if the time-dependence 
of a and e are taken into account. This integration 
was continued just until the IMBH/SMBH binary 
had coalesced. 



system is highly chaotic and the figure shows that 
in this very difficult case the new method is more 
accurate and much faster. 

Figure 2 shows a second numerical experiment 
that included all post-Newtonian terms up to or- 
der PN2.5. We considered the problem of peri- 
astron shift of a single star orbiting around the 
supermassive black hole (SMBH) at the center of 
the Milky Way. The particle masses were mi = 
3.5 X IO^Mq and m2 = lOM© = 2.85 x lO^^mi, 
and the Keplerian orbit had initial semi-major 
axis a = 0.01 pc, similar to that of the "S " 



stars ( Ghez et all 120051 : lEisenhauer et all boOSh . 



In units where G = mi = 1, and adopting 1 
mpc = 10~^ pc as the length unit, the speed of 
light is 77.19. Three different eccentricities were 
tried: (0.9,0.98,0.99). The same values of (a, 7) 
were used as in the first experiment. The integra- 
tions were continued for 10^ yr or ~ 20, 000 orbital 
periods. The figure plots the orientation of the 
Laplace-Runge-Lenz vector every 1000 orbits; the 
solid (red) lines show the periastron advance ex- 
pected based on the PN2.5 equations, which pre- 
dict a shift each period of 



6ttGMi 



a(l-e2)c2 2a2(l-e2)2c4 



(73) 



or (0.095, 0.45, 0.91) degrees for e = (0.9, 0.98, 0.99). 

As a third experiment, we investigated the ef- 
fect of an intermediate-mass black hole (IMBH) 
on the motion of a star orbiting around the Milky 
Way SMBH. The IMBH was given a mass of 
3500Mq = 10~^MsMBH, semi- major axis 0.1 mpc 
and eccentricity 0.9. The large eccentricity greatly 
reduced the gravitational-wave inspiral time, by 
a factor ~ 10'^ compared with a circular orbit; 
inspiral required ~ 3670 yr or 10^ initial or- 
bital periods. To this binary system was added a 
third component of mass IOMq on an orbit hav- 
ing a = 8 mpc « 1600 AU and e = 0.974 with 
respect to the center of mass of the SMBH/IMBH 
binary. These orbital elements are similar to those 
inferred for the star SO-16 (Ghez et al. 2005). The 
star and IMBH were coplanar with aligned angular 
momenta. The star completed ^ 10^ orbits dur- 
ing the time of IMBH inspiral. Figure [3] shows 
the evolution of the star's semi-major axis and 
its orbital orientation during this time, and the 
configuration-space trajectory is plotted in Fig- 
ure m Initially, the apoastron distance of the 
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IMBH is ~ 0.19 mpc, similar to the periastron 
distance of the star, ~ 0.21 mpc, so that close in- 
teractions are allowed. The primary influence was 
found to be on the star's semi-major axis; the ec- 
centricity of the star's orbit changed only slightly. 
The star's periastron advancement was found to 
be well predicted by equation ([73|) if the time 
dependence of a and e were taken into account 
(Fig. [3|). A second integration without the PN 
terms confirmed that the IMBH itself contributes 
only slightly to the star's precession rate for this 
configuration. 

Our final numerical experiment was a one- 
million-year integration of a seven-body system 
consisting of the MW SMBH; an IMBH of mass 
10~'^Msmbh; and the five, shor test-period S stars : 
SO-1, SO-2, SO-16, SO-19, SO-20 ijOhez et al.ll20(ih . 
Stellar masses were set to ISM© and the initial 
positions and velocities were determined at year 
2000 AD using the Keplerian orbital elements 
given in Table 3 of Ghez et al. (2005). The 
IMBH orbit was assigned an initial eccentricity 
of 0.9 and semi-major axis of 1 mpc, compared 
with 4 mpc ^ a < 25 mpc for the stars. Post- 
Newtonian terms were included, causing the orbit 
of the IMBH to precess rapidly and to fill the an- 
nulus 0.1 mpc ^ r < 1.9 mpc; for this choice of 
{a,e, Mimbh) the gravitational wave inspiral time 
is ~ 10^ yr, much greater than the length of the 
integration. Three of the included stars, SO-2, SO- 
16, SO-19, have periastron distances that intersect 
the IMBH's orbital annulus and so each of these 
stars was able to interact closely with the IMBH, 
although many orbital periods were required be- 
fore close encounters occurred. SO-19 was the first 
to achieve positive energy, at t w 147, 500 AD; 
before ejection the star's orbit evolved toward 
small a and e. SO-16 was the next to escape, at 
t « 254, 500AD; this star moved into a highly ec- 
centric orbit before being ejected. SO-2 remained 
bound to the SMBH/IMBH binary but its semi- 
major axis increased gradually to ~ 1 pc, roughly 
equal to the radius of influence of the SMBH (if 
the latter were embedded in the Galactic bulge), 
so in a practical sense it too may be considered to 
have escaped. Figure [5] shows the evolution of a 
and e for the three stars. Experiments like these 
could be used to constrain the mass and orbital 
parameters of a putative IMBH near the Galactic 
center. 



6. Concluding Remarks 

We have demonstrated the logH method (with 
leapfrog and Bulirsch-Stoer-extrapolation) is at 
least equally good, and for some systems better, 
than the KS-CHAIN method. The optional TTL 
(which we use only to provide some regularization 
for the occasional close encounters of very small 
bodies) seems to have some problems that may 
have to do with roundoff error when the func- 
tion uj first increases then decreases by a large 
amount in case of very close encounters. How- 
ever, we point out that this round-off problem is 
reduced almost to non-existence by including in 
the TTL-function f2 only the interactions between 
small bodies and using only a small factor for 
this quantity in the time transformation function. 
Since the TTL function derivative is explicitly in- 
tegrated it automatically reduces the stepsize in 
case of close approach of small bodies thus provid- 
ing more careful integration of such events even if 
the coefficient f3 is negligible. This is due to the 
properties (sensitivity) of the Bulirsch-Stoer ex- 
trapolation method that sees any tiny irregularity 
in the data and modifies the stepsize in response. 
This property of the BS-extrapolation method also 
helps if the leapfrog algorithm happens to evaluate 
the derivatives too close to collision of some pair 
of bodies: the step is rejected and re-evaluated 
with a different stepsize. This procedure normally 
avoids the repetition of a such a numerical acci- 
dent. Finally, our code gives the option of not 
using time transformation at all. Even in this al- 
ternative the present code is much better than the 
straightforward use of center-of-mass coordinates. 
The reason is that the chain coordinates reduce 
the roundoff significantly. This is one of the great 
advantages of the chain structure. 



DM was supported by grants AST-0420920 and 
AST-0437519 from the NSF and grant NNXG7AH15G 
from NASA. 
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X (mpc) 

Fig. 4. — Configuration-space trajectory of the 
star whose orbital elements are plotted in Fig. [3l 
The orbit remains in the X — Y plane; the semi- 
major axis is initially parallel to the X axis. The 
semi-major axis changes randomly due to pertur- 
bations from the IMBH/SMBH binary, but after 
the latter has shrunk appreciably, a remains nearly 
fixed and the only evolution is uniform precession 
due to the post-Newtonian terms. 
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t (yr) 

Fig. 5. — Evolution of the semi-major axes (top) 
and eccentricities (bottom) of stars SO-2, SO-16 
and SO-19 with respect to the Milky Way super- 
massive black hole in a 7-particle integration that 
included an intermediate-mass black hole and the 
five, shortest-period SO stars. Time zero corre- 
sponds to 2000 AD. Arrows in the top panel in- 
dicate when SO-19 and SO-16 are ejected; SO-2 re- 
mained formally bound to the SMBH/IMBH bi- 
nary but its semi-major axis gradually increased 
to ^ 1 pc. 
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